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Abstract 

We hereby introduce and study an algorithm able to search for initial con- 
ditions corresponding to orbits presenting forced oscillations terms only, 
namely to completely remove the free or proper oscillating part. 

This algorithm is based on the Numerical Analysis of the Fundamental 
Frequencies algorithm by J. Laskar, for the identification of the free and 
forced oscillations, the former being iteratively removed from the solution 
by carefully choosing the initial conditions. 

We proved the convergence of the algorithm under suitable assumptions, 
satisfied in the Hamiltonian framework whenever the dAlembert character- 
istic holds true. In this case, with polar canonical variables, we also proved 
that this algorithm converges quadratically. We provided a relevant appli- 
cation: the forced prey-predator problem. 
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1. Introduction 

A large number of low-dissipative problems can be conveniently modelled 
over a quite short time interval as if they were conservative systems, i.e. in 
neglecting the short term influence of the dissipative processes, and more- 
over to evolve close to a stable equilibrium because the oscillations around 
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it can be assumed to have been damped in the past. This is often the case 
in celestial mechanics, where for instance the resonant rotation of planetary 
satellites is assumed to be at an equilibrium, and thus only depending on an 
external torque, without any free oscillation. However, the problem is more 
generally applicable to a large class of quasiperiodically forced systems, for 
instance we will provide an application to a well-known prey-predator sys- 
tem subjected to an external periodic forcing, accounting for the seasonal 
changes. Let us observe that the method is far more general and the appli- 
cability domains goes well beyond the two examples hereby presented as it 
can be seen. 

More precisely we consider systems having stable equilibria and quasi- 
periodic orbits in their neighborhoods. We start from a system defined by 
an ordinary differential equation alike 

X{t) = f{X)+g{X,t), (1) 

where / : B C C" ^ C", 5 : B C C" x M ^ C" being an external per- 
turbation, and X £ M G C", B being an open subset of C" containing the 
equilibrium we are interested in, and we assume that every component Xi{t) 
of the solution, X{t), in a neighborhood of the equilibrium can be expressed 
as a convergent infinite sum of periodic trigonometric monomials, such as: 



= ^a«exp {if it) 



where a/ are complex amplitudes and // are real frequencies. We can split 
these frequencies into 2 groups: free, or proper frequencies of the "unper- 
turbed" system, i.e. when g = 0, and forced frequencies due to the external 
perturbation g. 

This decomposition holds for quasi-integrable systems in the Hamil- 
tonian framework thanks to the Kolmogorov-Arnold-Moser (KAM) theo- 
rem Esl . 22|. that ensures the existence of quasi-periodic invariant tori 



under suitable hypotheses: the non-degeneracy of the integrable part of the 
Hamiltonian function, the Diophantine condition on the frequencies and the 
smallness of the perturbation; let us observe that our setting is a simpli- 



fied version of the Poschel result 3l| where the normal frequencies, hereby 



called proper frequencies, do not depend on the torus frequencies, hereb 



named forced ones. On the other hand, the theorem by Nekhoroshev |25l.l26l| 
proves that the orbits of a slightly perturbed system stay close to the orbits 
of the unperturbed system over an exponentially long time with respect to 
the inverse of the perturbation, and thus ([2]) is a good approximation for 
any realistic physical times. 
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For example, in the context of the resonant rotation of the Moon around 
the Earth, we can spht the frequencies fi into two groups: the frequencies uj 
that are forced perturbations displacing the equilibrium, and the frequencies 
Ui that are free librations of the system around the equilibrium. The fre- 
quencies i^j are due to an external forcing (in our example, the gravitational 
perturbations of the Sun and the planets on the rotation of the Moon), while 
iOi depend on the intrinsic properties of the system (in our example, on the 
gravity field of the Moon). Slow dissipative processes acting since the origin 
of the system (4.5 Gyr in the present case) are expected to have damped 
the free librations (see e.g. [sO] for the rotation of Mercury, and 12] for the 
planetary satellites). 

The numerical computation of the wanted orbits, namely whose free li- 
brations have negligible amplitudes, i.e. the most realistic ones, requires 
the accurate knowledge of the initial conditions corresponding to the equi- 
librium. This task can be tricky whenever the system is subject to several 
different sinusoidal perturbations. In such a case, an analytical determina- 
tion of the equilibrium would either give a too approximate solution because 
related to a too simple system, or would be too difficult to use for the com- 
plete system to obtain the required accuracy. 

A straightforward way to overcome this difficulty is to use a numerical 
algorithm that iteratively corrects the initial conditions by identifying and 
then removing the free oscillations around the equilibrium. The idea of using 
the frequency analysis to refine the initial conditions is not completely new 
in the literature, in fact it has already been used by Locatelli & Giorgili 



2l[ for a computer-assisted proof of the KAM theorem. Here, in removing 



the free oscillations, we aim at reducing the numbers of dimensions of the 
solutions, i.e. we are looking for solutions lying on lower dimension tori with 
respect to the full dimension n+p. This has already been done by Noyelles 



et al. (see e.g. [Ul, |28|, |29|) and Robutel et al. [3j] for representing the 
resonant spin-orbit rotation 3:2 or 1:1 in the Solar System, by Couetdic et 
al. [3] in the framework of (exo) planetary systems, and by Delsate (3] for 
the dynamics of a spacecraft around the asteroid Vesta. The aim of our 
paper is to clearly state the algorithm, discuss its applicability domains and 
moreover prove under suitable hypotheses its convergence. 

The paper is organized as follows. Section [2] is devoted to the presenta- 
tion of the algorithm and the numerical tools needed. Then in Section [3] we 
will present a convergence proof of the method in the Hamiltonian frame- 
work, before addressing its extension to a more general framework in Sec- 
tion m In Section [5] an application of the algorithm will be presented, to 
a problem from mathematical biology. We finally review briefly other al- 
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gorithms in Section [6l before summing up and drawing our conclusions in 
Section [71 



2. The algorithm 

The goal of this section is to present an algorithm to search for the 
equilibrium of a periodically forced system. Let observe that its aim is not 
to prove the existence of the wanted orbit that should be assumed to exist 
a priori, but to determine the initial condition providing the wanted orbit. 
For the conditions of existence of such orbits, we refer the interested reader 
to e.g. Hi]. 

The algorithm requires hence the existence of quasi-periodic orbits around 
the equilibrium and thus it relies on a method able to reconstruct these 
quasi-periodic orbits by identifying the forced and the free terms in the fre- 
quency space. To accomplish this task, in the following we decided to use 
the NAFF alg orithm, Numerical Analysis of the Fundamental Frequencies, 



by Laskar [18|, [19|, shortly described in Appendix A. We used this algorithm 
because it has proved its efficiency, but any other algorithm allowing to 
detect accurately the frequencies could be used as well. 

Let us consider a system described by the differential equation 

^ = f{X) , (3) 

in a neighborhood of an either (quasi)periodic or static solution with fre- 
quencies {uJi)i<i<:m with < m < n written as cj G M", hereby named the 
equilibrium of the free motion, X being a n-dimensional vector. Let us then 
modify the system by adding a quasiperiodic forced term: 

J( = f{X) + g{X,t), (4) 

where g{X,t) is quasiperiodic, its frequencies being (i^j)i<j<p, also written 
as z7 E M^. So, the functions / and g are defined from an open set of C" and 
C" X M respectively, to C", and n and p are strictly positive integers. 

Let us assume that and Uj are rationally independent and thus secular 
terms are not allowed in the solution. Let us also denote the solution with 
a generic initial condition X by (j){t;X) : M x C" — )• C". Using NAFF we 
can identify the contribution of each frequency and thus obtain the decom- 
position 
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that can be rewritten by separating free from forced oscillating terms, as 
follows 



=: S{t;X) + L{t;X), (6) 

where the quantities S and L have been here defined by identifying left and 
right hand sides, they are functions defined from an open set of C" x M to 
C"". Their uniqueness comes from the assumption that the determination of 
the quasi-periodic series has a perfect accuracy. 

The assumption of the existence of a {2t: /vj)-^^.^^ quasiperiodic orbit, 
i.e. lying on a p-torus, namely composed only by forced oscillating terms, 
translates into the existence of an isolated initial condition such that the 
solution with initial condition X^o is {'^'^ /^j)i<^j<^p quasiperiodic, namely: 



{t] ^oo^ — S{t] Xoo) ) 



or equivalently L(t;Xoo) = 0. X^o is unknown and our algorithm is a way 
to determine it. 

This algorithm can be formulated as follows: 

1. Take Xq sufficiently close to Xoo] 

2. Integrate the ODE (jl]), thus determine the solution cf (^t; Xq^ and then, 
using NAFF, get the decomposition (i)(^t;Xo^ = S{t;Xo) + L{t]Xo); 

3. Define Xi = S{0;Xq) and iterate point 2 using Xi instead of Xq] 

4. In this way the algorithm will produce a sequence Xn, iteratively de- 
fined by 

Xn+i = S{0-Xn), (7) 

such that, when n — )• cxd, Xn — t- Xqo, or equivalently L{t; Xn) — s- 0. 

Remark 1. Once numerically implemented, we can define at least two ways 
to define a stopping criterion for this algorithm: we can consider the con- 
vergence as reached when the largest amplitude associated with the free os- 
cillations is too small to be detected, or when it is small enough to not 
significantly alter the determination of the forced oscillations. 
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We here make no hypothesis on the dimensions n and p. This algorithm 
has aheady been used in the following contexts: 



• Noyelles (2009) [2^]: n = 3, p = 13 

• Dufey et al. (2009) n = 2,p = 5 

• Couetdic et al. (2010) 0: n = 1, p = 3 

• Noyelles et al. (2010) n = 4, p = 5 



Robutel et al. (2011) [3|]: n = 1, p = 6 



• Delsate (2011) n = 1, p = 3 

Its convergence under some assumptions will be proved in the next two 
sections. 

3. Convergence of the algorithm in the Hamiltonian framework 

Let us consider a (n + p)-DOF Hamiltonian H{yi,Xi,Aj,\j), Ui {1 < i < 
n) and Aj (1 < j < p) being the momenta and Xi and \j the variables. 
We assume that the forcing is due to Aj,Xj and that they are respectively 
actions and angles variables. Moreover, we hypothesize that the Hamiltonian 
system can be locally described by a forced perturbed harmonic oscillator: 

n{Ui,Ui,Aj,Xj) = '^UiUi + eni{Ui,Ui,Aj,Xj) , (8) 

i 

where is the frequency of the free oscillations and eTii is a small pertur- 
bation. This is a degenerate setting where existence of invariant tori can be 
studied thanks to Birkhoff theory (see e.g. [2^). The Hamiltonian T-L has 
been obtained after canonical transformations among which an untangling 



one that removes the second-order cross terms alike yji/k with j ^ k 14l |. 
and the classical canonical polar transformation: 



Xi = y/2UtZi sin Ui, 

Vi = y/2Ui/Zi cos Ui, (9) 

where Zi is a constant term, suitably chosen so that the first order terms in 
\/TJi disappear. The transformation ([9]) can be seen as the expression of the 
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variables Xi and yi after averaging with respect to the forced perturbation, 
i.e. 



\xi{Ui,Ui,-,-) = ^/2UiZiS\nui, 
\yi{Ui,Ui,-,-) = ^JlVilZiCosui . 

There are several perturbative methods, based on the hypotheses that 
the perturbation is small and far enough from resonances with the proper 
frequencies Wj, that allow to derive the complete expression of Xj and y,, i.e. 
including their dependencies on Aj and A , . We hereby propose to use the 
method of the Lie transforms (see or 10] for detailed explanations). In 



general, the formal series do not converge, nevertheless, such perturbation 
techniques can be justified through Poincare's theory of asymptotic approxi- 
mations [s^]- It is the difference between the convergence seen by geometers, 
on a finite-size interval, and the one seen by astronomers, who see the series 
as asymptotic expansions around 0. 

Under the assumption of an analytic Hamiltonian Ji, Henrard fis'] has 
proved that at the A^-th step of the Lie transforms algorithm, the functions 
Xq^^ and i/g^AT with 1 < g < n, follow the d'Alembert characteristic for 
{Ui^N,Ui,N), i-e. they can be written as 



IM / 

ki,hj& ki m>l 



X exp [i[ ^ kiUi^N + ^ Aj 

(10) 

and 

ki,hj£Z ki m>l 

X exp (^(y^feitti,Af + Aj) j. 



(11) 

where ak^^hp Pk„hp 1m,i,ki,hj and S^m^kijhj 'ire complex constants, and rrii, 
ki and hj are integers. The right-hand members are assumed to converge 
absolutely for Ui^N < U' for some positive U' . 
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So we get the time evolution of the above functions given by: 



Xg,N 



ki,hj€Z ki m>l 

X exp (i[ E kiUi^N{t) + E hjXjit) 



(12) 



ki,hjEZ ki m>l 

X exp E hui,N{t) + E ^J-^J (*))) ■ 



(13) 



Assuming that the truncation order M is large enough such that the 
functions f/j,jv(t) ~ U'^ ^ is almost constant, as they would be in the limit 
M — )■ oo, we can obtain the decomposition of (t), and similarly for (t), 
into forced and free oscillating terms as follows 



ki,hj& m>l hj 
^ ' 

Sq^N 

ki^OjhjEZ ki m>l ki hj 



(14) 



For any A'' the Sq^M terms, components of S'at behave as 



S<l,N ~ ^ + E Bi,j\/KN\/U'j,N + ■■■ 

while the Lg^N terms, components of Lg^N, behave as 

i 

, where " . . . " means " higher order terms" . 



Let us observe that setting t = in the Eq. [T3] we get from the very 
definition of our algorithm: 

Xq,N+l = Sq^NiO) , (15) 

that should equals the left hand side of O written with + 1, that is: 

x,,^+i ~ A + ^ B,j^ful^^fu[^ + ... (16) 

From the EqJ15landll6lwe get 

^Ul^^, = 0{U,,n). (17) 

If we assume the different components Ui^n of the vector Un at the iter- 
ation to be of the same order of magnitude, then we have a quadratic 
convergence of the algorithm. 



4. The convergence in a more general framework 

In a general framework, the convergence cannot be checked a priori. We 
here investigate a condition leading to this convergence, before discussing 
its relevance. 

4-1. The hypothesis 

The convergence of the algorithm can be proved by showing the conver- 
gence of the sequence (XA:)fceN defined by ([7]). 

If we have — s- P for some P, and P being of dimension n, then 
by rewriting ^ as follows 

Xk+i = 5(0; Xk) = 0(0; X^) - L(0; X^) = X^ - 1(0; Xk) , (18) 

we straightforwardly get by continuity of L and the convergence hypothesis 
oiXk 



= lim X„+i - X„ = - lim L(0; X„) = -L(0; P) , 
then by uniqueness of the 27r/i/j-quasiperiodic orbit we can conclude that 

P = Xoo- 

We provide a proof of the convergence of the algorithm assuming the 
following hypothesis: the Jacobian matrices E and A of the functions S and 
L do satisfy 



lim 



-1 



A{X) S(A:) = 0. (19) 
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4.2. Proof 

From the definition ([7]) of X^+i and ttie very definition of X^+i as initial 
condition of the orbit we get: 

5(0; Xk) = Xk+i = 0(0; Xk+i) = 5(0; Xk+i) + HO; ^fc+i) • (20) 

This relation defines implicitly the map that associates to Xk the next iter- 
ation Xk+i- Let us denote it for short by F, namely Xk+i = F{Xk). 

Assuming the existence of a 271 /uj quasiperiodic orbit corresponding to 
the initial condition X^o is equivalent to assume that F admits X^o as a 
fixed point, i.e. F ^Xqo^ = ^oo- Then, to prove that X^ converges to X^o 
we have to prove that this fixed point is indeed an attractor, namely that 
every eigenvalue of the Jacobian matrix $ of F has a modulus lower than 1. 

Rewriting (|20p as 

5(0; F{Xk)) + L(0; F{Xk)) = 5(0; Xk) , (21) 

and then differentiate every component Si of it with respect to every com- 
ponent Xi of Xk'. 

dS^^^P_ f dF{Xk)j dSi ^ dFiXk), dLi\ 
dxi \ dxi dxj dxi dxj I 



In calling S the Jacobian matrix of 5, we straightforwardly get from (j22p 
S(Xfc) = S [F{Xk)) HXk) + A {F{Xk)) HXk), (23) 



what yields 



$(Xfe) = (s(F(Xfc)) +a(f(X,,))) 's(Xfe) (24) 



j:f^F{Xk))+A[F{Xk))) (25) 



xA (F(Xfc) [a (F(Xfc)] ' S(Xfc) (26) 
From the hypothesis (|19|) we have 

Jiin $(Xfe) = 0, (27) 

SO the eigenvalues of the matrix <I>(Xoo) have a modulus lower than 1. 
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4-3. Discussion 

The caveat of this proof is that the hypothesis (jl9p usuahy cannot be 
checked a priori, so for a given dynamical system studied close to a stable 
equilibrium, we cannot a priori know whether our algorithm will converge 
or not. Since this algorithm aims at determining the solution P = X^o, it is 
a priori unknown as well. 

Anyway, we can prove that the hypothesis ()19p is verified in the Hamil- 
tonian framework, since it is in fact a weakening of the d'Alembert charac- 
teristic. For sake of clarity, we here reduce to 1-degree of freedom systems, 
but the principle is the same for multi-dimensional systems. In this context, 
the hypothesis (fT9|) reads 

lim|#4=0^ (28) 
Thank to the d'Alembert characteristic we have: 



S'(0; x) ~ + — a^ool + • • • and L{0; x) ~ by^\x — Xoo\ + . . . . 

Thus dxS{0;x) is finite while dxL{0;x) diverges, as y^|x — Xoo|, and condi- 
tion (fT9|) is satisfied. 

More generally in the case 

S{0; x) ^ Xoo + a{x — Xao)" + ■ ■ ■ and L(0;x)~ b{x-x^f + ... , (29) 



for some positive a and /?, such that a > (3 satisfies the condition (jl9p . one 
can provide also the rate of convergence of the algorithm. In fact, let us 
define On = Xn — Xoo, then using ([29]) . Eq. (pOj) rewrites: 

that can be approximately solved for ^n+i, for instance using the Lagrange 
inversion formula [5], to get 



and thus providing a convergence rate a/ j3. In the above mentioned Hamil- 
tonian framework this results into a quadratic convergence rate. But to have 
the convergence of the algorithm, we just need (3 > a. 
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5. Application of the algorithm to a problem in mathematical bi- 
ology 

In this section we want to investigate the appUcability of the algorithm 
beyond the Hamiltonian framework. We thus chose an example from math- 
ematical biology, a prey-predator system, derived from the Lotka-Volterra 
equations, periodically forced by a sinusoidal term, accounting for the sea- 
sonal changes. 

5.1. Forced prey-predator systems 

This model presented in Blom et al. [3] aims at representing the densities 
of preys xi and predators X2 as a function of time; beside the standard inter- 
actions among predator and prey and the logistic growth rate for the prey, 
we also take into account an external periodic forcing term, that can rep- 
resent the one-year-periodic density of food availability due to the seasonal 
alternation. 

The equations ruling the densities of the populations read: 

axi{l + ^ cos{2'Kt) — X2 — r]Xi) , (30) 

/?X2(-1 + Xi), (31) 

where a, /3, 7, 5 and 77 are non negative constants. When 7 = (autonomous 
system) and 77 = 0, one can straightforwardly prove the existence of T = 
27r/\/a/3 periodic solutions, around the non trivial equilibrium {xi,X2) = 
(1, 1). If 77 > the non trivial equilibrium moves to (xi, X2) = (1,1 — f]) 
moreover these oscillations are damped; for 7 > the system is subjected 
to 27r-periodic forced oscillations (see Fig. [T] left panels). 

Using the algorithm, we can determine the initial conditions correspond- 
ing to the solution with zero amplitude (see Fig. [T]right panels, black curves) 
associated with the proper period T, without using any damping (i.e. with 
r? = 0). 

5.2. Numerical analysis 

We performed a numerical analysis of the above forced predator-prey 
model using the following parameters values: a = 4.539, (3 = 1.068, 7 = 0.25 
and r/ = 0. We adopted the variable step size Bulirsh-Stoer algorithm 0, 113] 
to numerically integrate the differential equations (f30|) and (|3T]l . The proper 
frequency, for the case 7 = 0, is given by \Jaf3 = 2.201 74. Results are 
reported in Table [H 



dxi 
IT 

dX2 

IT 
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Ti=0.025 Ti=0 




t t 



Figure 1: Numerical integration of the system of Blom et al. for a = 4.539, = 1.068 and 
7 = 0.25, the initial conditions being xi = 1 and X2 ~ 1 — rj (in black), and the period 
T = 2.853 74. We can sec a damping of the oscillations for r/ = 0.025. On the right, the 
grey curves result from numerical integrations of the equations of the system with initial 
conditions obtained after 2 iterations (n = 2). The resulting behavior is very similar to 
the one obtained after damping. 

Table 1: Application of the algorithm to the forced prey-predator system, with a = 4.539, 



P = 1.068, 7 = 0.25 and ?? = 0. 
n I.e. 


#freq. Rank 


Ampl. io* 


.7:1 1.000 000 000 000 000 
.V2 1.000 000 000 000 000 


50 2 
50 1 


3.831 163 X 10-^ 2.206 634 
1.834 280 X 10"- 2.20GG34 


xi 0.989166 714 745100 
X2 0.965 514 795157481 


22 4 
21 3 


3.573 335 x 10--'^ 2.207483 
1.729 063 x 10-5 2.207483 


XI 0.989186 585 234344 2971 
X2 0.965 545142 090197 513 7 


26 6 

26 6 


4.508 632 X IQ-y 2.207483 
2.181634x 10-9 2.207483 


xi 0.989186 576 347806470 6 
X2 0.965 545142191326 750 5 


14 11 
14 11 


6.524090 x 10-" 2.207483 
3.156 872 x 10-1^ 2.207483 


xi 0.989186 576 347 806 470 2 
X2 0.965 545142191326 7504 


14 11 
14 11 


6.503088 x 10-^" 2.207483 
3.146 710 x 10-" 2.207483 
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We see variations of the detected proper frequency a;*, before converging 
to 2.2074, that is close to the predicted value of 2.201 74 (with 7 = 0). 

In right panels of Figure [U we plot (grey curves) the evolution of the 
variables xi and X2 for the initial condition corresponding to the iteration 
n = 2 (Table [TJ. We notice that, with these initial conditions, we removed 
the free libration parts of the signal. The resulting behavior is very similar 
to the damping case (rj ^ 0) but with the damping the equilibrium is shifted 
from (1,1) to (1,1 - r]). 



6. Alternative algorithms 

In this section we briefly review other existing algorithms having goals 
similar to the one we present here. We do not intend to make a rigorous and 
exhaustive comparison between them and the presented algorithm because 
we think that these studies go beyond the scope of this paper. We thus 
deserve this for forthcoming papers. We here focus on 2 algorithms: the 



iterative method by Rodionov et al. ( 38|, |37|, |35|, l36|] ) , and the fit of libration 
centers by Bois & Rambaux [1]. 

6.1. The iterative method (Rodionov et al. 2006-2011) 

This method has been elaborated in the framework of galactic dynam- 
ics, to find the initial conditions of N-body simulations corresponding to 
equilibrium states. These equilibrium states are characterized by various 
parameters that describe for instance the mass density of a galaxy, or the 
radial velocity dispersion. This algorithm consists in choosing a priori suit- 
able initial conditions, letting the system evolve in propagating its orbit, 
possibly under some constraints, and use the final state of the system to 
build new initial conditions. 

The authors claim that this method can be adapted to any arbitrary 
dynamical system, and we agree with this opinion, at least for the ones 
modeling physical phenomena. A rough adaptation of this algorithm in the 
context of the rotation of celestial bodies in spin-orbit resonance would con- 
sist in including a dissipative effect in the equations of the problem, and 
letting the system evolve on a long enough timescale (Peale [30] estimates 
this timescale of the order of 10^ years for the rotational dynamics of Mer- 
cury) for the damping to act efficiently, and then to use the computed spin 
angle and spin velocity to start a new simulation. A way to bypass the 
problem of the long damping time would be to artificially accelerate the 
damping process. Nevertheless, the dissipation changes the equilibrium po- 



sition of the system (see e.g. [SSy). A refinement of this process could be 
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to use a fast dissipation at the first iteration, and after a slower one. Using 
the presented algorithm consists in making the assumption that the effect 
of the dissipation on the location of the equilibrium is negligible. In fact, 
this effect has not been observed yet except for the Moon; detecting it for 
Mercury with the space missions MESSENGER and Bepi-Colombo is chal- 
lenging. So, this assumption can be considered as acceptable. An advantage 
of the iterative method is that it does not assume the existence of a 2TT/vj 
quasiperiodic orbit. 

In the case of the prey-predator problem (SeclJ]), ths algorithm requires 
to set the parameter r/ to 0, while the iterative method would use its strictly 
positive value. As we have already seen, it shifts the mean equilibrium from 
(1,1) to (1,1-r?). 

6.2. Fit of libration centers (Bois & Ramhaux 2007) 

This algorithm has been elaborated in the context of a numerical compu- 
tation of Mercury's equilibrium obliquity, and to the best of our knowledge 
has neither been used in any other study. This problem is a 2-dimensional 
spin-orbit problem of a rigid body, with the notable difficulty that it con- 
tains a very long period due to the regressional motion of Mercury's orbital 
ascending node, whose period is of the order of 250, 000 years. An efficient 
numerical identification of this sinusoidal term would require to propagate 
the equations of the system over at least 5 x 10^ years (i.e. two periods), 
while the space missions cover a period of about 2 years, and the orbital 
ephemerides available cover a few thousands of years. So, trying to iden- 
tify the period of regression of the node seems not to be the right way to 
do. Anyway, the use of unoptimized initial conditions induces 1,000-y free 
oscillations whose amplitude is expected to be null in the real system. So, 
Bois & Rambaux propose to fit a sinusoid to the orbits obtained after nu- 
merical simulations (in which the dissipation is neglected) to identify the 
free oscillations, and to remove them from the initial conditions to perform 
a new simulation. This way of identifying free oscillations without making a 
complete frequency analysis of the system does not require the final solution 
to be 2iT/vj quasiperiodic (in fact, the 250,000-y forced oscillation is seen 
as the slope over this timescale). So, this is just a partial decomposition of 
the signal, a complete one, when possible, is of course more accurate. 

7. Conclusion 

We hereby presented an algorithm able to determine the initial condi- 
tions corresponding to forced oscillating equilibria, namely removing free 
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proper librations. In our restricted group of local collaborators, we are used 
to name it NAFFO, for Numerical Algorithm For Forced Oscillations, since 
this recalls its proximity to NAFF. But this is just for internal use, we do 
not claim for the fatherhood of this algorithm and every user is of course free 
to give it the nickname he prefers. Given an initial condition, this algorithm 
iteratively produces better and better approximated initial conditions, by 
integrating the system, identifying the free and forced frequency terms and 
eventually removing the former ones. In the present paper this step has 
been performed using the NAFF method by Laskar. 

Under suitable conditions, we proved the convergence of the algorithm, 
under the assumptions of exact numerical integration and frequency iden- 
tification. We shown that in the Hamiltonian framework the required hy- 
pothesis are satisfied, whenever the d'Alembert characteristic holds true. 

We provided one relevant application to benchmark the proposed algo- 
rithm, in mathematical biology. This supplements the applications already 
present in the literature. 
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Appendix A. The NAFF algorithm 

Since the algorithm we present requires an accurate quasi-periodic rep- 
resentation of the orbit, we hereby shortly present the method we used: the 
NAFF (Numerical Analysis of the Fundamental Frequencies) algorithm due 



to J. Laskar (see e.g. [18], [191] where its convergence has been proved under 
suitable Diophantine hypothesis for the frequency vector). The good accu- 
racy properties of the NAFF, allows to use it for instance to characterize the 
chaotic diffusion in dynamical systems El in particles accelerator 24 1 or 
in celestial mechanics [27j, |20(] . 



The algorithm aims at numerically identifying the coefficients ai and ui 
of the quasi-periodic complex signal x(t) known over a finite, but large, time 



In this framework, it is also known as FMA for Frequency Map Analysis. 
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span [—T;T], of the form ([2]), thus providing an approximation 

TV 

x{t)^^a:exp{^f;t) , (A.l) 
1=1 

where f*, respectively a*, are the numerically (as emphasized by the bullet) 
determined real frequencies, respectively complex coefficients, i.e. ampli- 
tudes. If the signal x{t) is real, its frequency spectrum is symmetric and the 
complex amplitudes associated with the frequencies i^' and — f* are complex 
conjugates. 

The frequencies and associated amplitudes are found with an iterative 
scheme. To determine the first frequency v*, one searches for the maximum 
of the amplitude of 

= {x{t),exp{iiyt)) , (A.2) 
where the scalar product {f{t),g{t)) is defined by 

{fit),g{t)) = ^ f_J{t)W)x{t)dt, (A.3) 

where g{t) is the complex conjugate of g{t) and where x(^) is a weight 
function, i.e. a positive function verifying 

^ j\{t)dt = l. (A.4) 

Laskar advises to use 

x(t) = -^(l + cos(7rt)j , (A.5) 

where p is a positive integer. In practice, the algorithm is the most efficient 
with p = 1 or p = 2. We used p = 2 in the numerical applications presented 
later in the paper, because it yields good results. 

Once the first periodic term exp(zz^*t) is found, its complex amplitude 
a\ is obtained by orthogonal projection, and the process is started again on 
the remainder fi{t) = f{t) — a* exp(zi/*t). The algorithm stops when two 
detected frequencies are too close to each other, because this alters their 
determinations, or when the number of detected terms reaches a maximum 
set by the user. When the difference between two frequencies is larger than 
twice the frequency associated with the length of the total time interval, the 
determination of each fundamental frequency is not perturbed by the other 
ones. Once the frequencies have been determined, it is also possible to refine 
their determination iteratively, numerically @] or analytically ||39l] . 
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